Quantitative analysis of numerical estimates for the permeability of porous media 

from lattice-Boltzmann simulations 



O 
(N 

-(— > 
o 

O 

m 



q 

O 

Oh. 



> 

o 
m 
in 

o 
o 



^ 



Ariel Narvaez,^'^ Thomas Zauner,'^ Frank Raiscliel,^ Rudolf Hilfer,^' ■^ and Jens Harting^-^ 

' Department of Applied Physics, TU Eindhoven, 

P.O. Box 513, NL-5600MB Emdhoven, The Netherlands 

"^ Institute for Computational Physics, University of Stuttgart, 

Pfaffenwaldring 27, D-70569 Stuttgart, Germany 

■^ Institute for Physics, University of Mainz, D-55099 Mainz, Germany 

(Dated: October 14, 2010) 

During the last decade, lattice-Boltzmann (LB) simulations have been improved to become an 
efficient tool for determining the permeability of porous media samples. However, well known 
improvements of the original algorithm are often not implemented. These include for example 
multirelaxation time schemes or improved boundary conditions, as well as different possibilities to 
impose a pressure gradient. This paper shows that a significant difference of the calculated perme- 
abilities can be found unless one uses a carefully selected setup. We present a detailed discussion of 
possible simulation setups and quantitative studies of the infiuence of simulation parameters. We 
illustrate our results by applying the algorithm to a Fontainebleau sandstone and by comparing our 
benchmark studies to other numerical permeability measurements in the literature. 
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I. INTRODUCTION 

The accurate numerical simulation of fluid flow in 
porous media is important in many applications rang- 
ing from hydrocarbon production and groundwater flow 
to catalysis and the gas diffusion layers in fuel cells [1]. 
Examples include the behavior of liquid oil and gas in 
porous rock [2], permeation of liquid in fibrous sheets 
such as paper [3] , determining flow in underground reser- 
voirs and the propagation of chemical contaminants in 
the vadose zone [4, 5], assessing the effectiveness of leach- 
ing processes [6] and optimizing filtration and sedimenta- 
tion operations [7] . An important and experimentally de- 
terminable property of porous media is the permeability, 
which is highly sensitive to the underlying microstruc- 
ture. Comparison of experimental data to numerically 
obtained permeabilities can improve the understanding 
of the influence of different microstructurcs and assist in 
the characterization of the material. 

Before the 1990's the computational power available 
was very limited restricting all simulations either to small 
length scales or low resolution of the microstructure. 
Shortly after its introduction lattice-Boltzmann (LB) 
simulations became popular [8-10] as an alternative to a 
direct numerical solution of the Stokes equation [11, 12] 
for simulating fluid flow in complex geometries. Histori- 
cally, the LB method was developed from the lattice gas 
automata [10, 13]. In contrast to its predecessor, in the 
LB method the number of particles in each lattice direc- 
tion is replaced with the ensemble average of the single 
particle distribution function, and the discrete collision 
rule is replaced by a linear collision operator. 

In the LB method all computations involve local vari- 
ables so that it can be parallelized easily [12]. With 
the advent of more powerful computers it became possi- 
ble to perform detailed simulations of flow in artificially 



generated geometries [3] , tomographic reconstructions of 
sandstone samples [8, 12, 14-16], or fibrous sheets of pa- 
per [17]. 

The accuracy of LB simulations of flow in porous media 
depends on several conditions. These include the resolu- 
tion of the discretization of the porous medium, proper 
boundary conditions to drive the flow and to implement 
the solid structure or the choice of the collision kernel. 
Even though advanced boundary conditions, discretiza- 
tion methods, as well as higher order LB kernels have 
been developed and are common in the literature, it is 
surprising to the authors that they only found limited 
applications so far. In particular for commercial appli- 
cations a three-dimensional implementation with 19 dis- 
crete velocities and a single relaxation time linearized 
collision operator is still the de-facto standard to calcu- 
late stationary velocity fields and absolute permeabilities 
for porous media [18] . Here, the flow is usually driven by 
a uniform body force to implement a pressure gradient 
and solid surfaces are generated by simple bounce back 
boundary conditions. 

The present work is motivated by the question whether 
permeabilities calculated by this standard LB approach 
can be considered to be accurate. In particular, it is im- 
portant to understand where the limits of this method 
are and how the accuracy can be increased. We quantify 
the impact of details of the implementation by studying 
3D Poiseuille flow in pipes of different shape and reso- 
lution and comparing the simulation results to analyt- 
ical solutions. This allows to demonstrate how simple 
improvements of the simulation paradigm can lead to a 
substantial reduction of the error in the measured perme- 
abilities. These include a suitable choice of the relaxation 
parameter r and the application of the multirelaxation 
time method in order to ascertain a minimal unphysical 
influence of the fluid viscosity on the permeability. Fur- 



thcr, a correct implementation of the body force to drive 
the flow together with suitable in- and outflow bound- 
aries is mandatory to avoid artifacts in the steady state 
velocity field. Finally, the small compressibility of the 
LB fluid requires a proper determination of the pressure 
gradient in the system. If these details are taken care of, 
it is shown that the LB method is well suitable for accu- 
rate permeability calculations of stochastic porous media 
by applying it to discretized micro computer-tomography 
(/i-CT) data of a Fontainebleau sandstone. 



II. SIMULATION METHOD 



The Boltzmann equation 



|./(x,c,0 



c.V/(x,c,t)=l](/(x,c,i)) (1) 



describes the evolution of the single particle probabil- 
ity density /(x, c,i), where x e IR'^ is the position vec- 
tor, c £ R"^ is the velocity vector, i € ]R is the time, 
and i7(/(x, c,t)) is the collision operator. While dis- 
cretizations on unstructured grids exists [19, 20], they 
are not widely used and typically the position x is dis- 
cretized on a structured cubic lattice, with lattice con- 
stant Ax. The time is discretized using a time step At 
and the velocities are discretized into a finite set of vec- 
tors Ci with i = 1, . . . ,N, called lattice velocities, where 
the finite integer N varies between implementations. In 
this work we exclusively use the so-called D3Q19 lat- 
tice, where iV = 19 velocities are used in a three dimen- 
sional domain [21]. A cubic lattice with basis ej, € M,^, 
/c = 1, 2, 3 is embedded into M,^ using the coordinate func- 
tion g : N^ !->• R^ to map the lattice nodes .£ G N^ to 
position vectors g(£) £ R'^. The computational domain 
is a rectangular parallelepiped denoted as 



C^{e£¥i^ ■.l<lk<Lk;k^l, 2, 3}, 



(2) 



where L^ £ N'^ are its dimensionless side-lengths. See 
Fig. 1 for a visualization. Physical quantities w such 
as pressure or density on the lattice are abbreviated 
as w{£.) = w{g{£.)). We introduce the vector notation 
f{£,t) = (/i(^,t), . . . , /Af(£,t)), where the components 
are the probabilities calculated as 

M£,t)^ f f /(x,c,i)dcdx. (3) 

JW{i) JB(i) 

Here, W(^) C R'^ is the finite volume associated with the 
point g(£) and B(i) C R'^ is the volume in velocity space 
given by lattice velocity Cj. The macroscopic density 
p{£,t) and velocity v(£,t) are obtained from fi{i,t) as 

N 

p{£,t) = p°Y,M£,t), (4) 

N N 

v{£,t) = ^ /.,(£, i)cj^/,(£,i), (5) 



where p° is a reference density. The pressure is given by 
p(£,t)=c,2p(£,t), (6) 

with the speed of sound [10, 13] 



1 fAx 



(7) 



Discretization of Eq. (1) provides the basic system of 
difference equations in the LB method 

fii£ + A£i,t + At)-fi{£,t)^Atni{£,t), (8) 

with A£i = CiAt/Ax and the initial condition fi{£,0) = 
1/N (for t ~ 0). The generally nonlinear collision opera- 
tor is approximated using the linearization 



ni{£,t) 



N 

E 



Si,{f,{£,t)~m£,t)), 



(9) 



probability function 
with a iV X iV 



around a local equilibrium 
r^{£,t)^ifP{£,t),...,f^^i£,t)), 
collision matrix S [22, 23]. 

The simplest approach to define the collision matrix 
uses a single relaxation time with time constant r. 



s.. 



^] 



1 . 



(10) 



where 5ij is the Kronecker delta. This single relaxation 
time (LB-BGK) scheme is named after the original work 
of Bhatnagar, Gross and Krook [24, 25]. Within the 
LB-BGK method, i°'-^{£,t) is approximated by a second 
order Taylor expansion of the Maxwell distribution [26] , 



.r(^,^)= 



p^c. 



1- 



2Cs4 



2cs 



(11) 



If external forces are absent, the equilibrium velocity is 
defined as v*(£,i) = 'v{£,t) from Eq. (5). As explained 
further below, v* (£, t) and v(.£, t) may differ from Eq. (5) 
if an external acceleration is present. The numbers uj^. 
arc called lattice weights and differ with lattice type, 
number of space dimensions and number of discrete ve- 
locities N . Sec [10] for a comprehensive overview on dif- 
ferent lattices. 

An alternative approach to specify the collision matrix 
is the multirelaxation time (MRT) method. Here, a linear 
transformation M is chosen such that the moments 



TV 



,(£,t) = ^A/i^/^(£,t) 



(12) 



represent hydrodynamic modes of the problem. We 
use the definitions given in [27], where mi{£,t) is the 
density defined in Eq. (4), m2{£,t) represents the en- 
ergy, mi{£,t) with i = 4,6,8 the momentum flux 
and 'mi{£,t), with i ~ 10,12,14,15,16 are components 
of the symmetric traceless stress tensor. Introducing 



the moment vector m(£,i) ~ {■mi{£,t), . . . ,mN{£,t)), 
Ct{£,t) = {ili{£,t), . . . ,^N{i,t)), a diagonal matrix 



Si 



h Si 



and the equiUbrium moment vector 



m^^il, t) = {mWl, t),..., rn^{l, t)), we obtain 



nie,t) 



M 



S-(m(£,t)-m°'i(£,t)). 



(13) 



During the cohision step the density and the momen- 
tum flux are conserved so that m'l^{£.,t) = mi{£.,t) and 
mi{£,t) = m'^^{£,t) with i = 2,4,6. The non-conserved 
equilibrium moments m1^{£^ i), i 7^ 1, 2, 4, 6, are assumed 
to be functions of these conserved moments and explicitly 
given e.g. in [27]. The diagonal element r, — l/sj in the 
collision matrix is the relaxation time moment rrii {£, t) . 
One has si = S4 = sg = ss = Oi because the corre- 
sponding moments are conserved, S2 = l/'T'buik describes 
the relaxation of the energy and sio — S12 ~ S14 = S15 = 
S16 = l/r the relaxation of the stress tensor components. 
The remaining diagonal elements of S are chosen as 

S = diag(0, 1/Tb,ik, 1.4, 0, 1.2, 0, 1.2, 0, 1.2, l/r, 

1.4, l/r, 1.4, l/r, l/r, l/r, 1.98, 1.98, 1.98), (14) 

to optimize the algorithm performance [27, 28]. Because 
two parameters r and Tbuik remain free, the multirelax- 
ation time method reduces to a "two relaxation time" 
(TRT) method. An alternative TRT implementation can 
be found in [29, 30]. 

To apply the LB method to viscous flow in porous me- 
dia it is necessary to establish its relations with hydro- 
dynamics. The Chapman-Enskog procedure shows that 
density, velocity and pressure fulfill the Navier-Stokes 
equations without external forces, with a kinematic vis- 
cosity [26, 31-34] 



v{t, At) = Cs^At 



T 

At 



Combining Eq. (15) and Eq. (7) gives 



T 

At 



VSi 



vAt 

'XAx^ 



(15) 



(16) 



Because v > Q^ Ax > 0, and At > 0, it follows that 
r/At > 1/2. 

A typical value for the pore diameter in sandstone 
is a ~ 10"^ m, and for water the kinematic viscos- 



ity and speed of sound are v 



10- 



1^ s ^ and Cs 



lO^ms"^, respectively. With typical velocities of order 
V ~ lO^^ms^^ the Reynolds number is Re = va/v w 
10"'^. Discretizing with Aa; = 10~^ m gives then r/At = 
0.5017. Because for r/At w 1/2 the LB method is known 
to be unstable, a direct simulation of water flow in porous 
media with these parameters is not feasible. To overcome 
this impasse, one might impose r/At = 1 and simultane- 
ously fix i> and Cs as fluid parameters. The discretization 
then is At w 10""'^^ s and Ax ss 10~^m. Again, a sim- 
ulation with these parameters is not possible because a 
typical pore with diameter a ~ 10~^ m would have to 



be represented by 10^ nodes, exceeding realistic memory 
capacities. Another way to circumvent these problems 
is to appeal to hydrodynamic similarity for stationary 
flows. The simulations in this paper are performed with 
fluid parameters that represent a pseudofluid with the 
same viscosity as water, but Cg — Ims"-'^ as the speed 
of sound. The discretization then is Ax = 10~^ m and 
At = 10~^ s. A pore of diameter a is then represented by 
10 nodes and a cubic sample with side-length 10"'^ m re- 
quires 1000^ nodes, a manageable system size on parallel 
computers. An external force, as discussed next, drives 
the flow such that the velocities are of order 10"'^ ms~^. 
The Mach and Reynolds numbers in the simulations are 
Ma « 10"'^ and Re w 10~^, characterizing a laminar 
subsonic flow. As long as Ma <^ 1 and hydrodynamic 
similarity remains valid, we do not expect that the pa- 
rameters of the pseudofluid will change the permeability 
estimate. 

An external acceleration h{£, t) acting on the fluid is 
implemented by adding two modiflcations. First, a forc- 
ing term written as a power series in the velocity [23] 



ipi{£,t) = At — ^|/io- 



P" 



2cs 



(17) 
is added to the right hand side of Eq. (8). Second, Eq. (5) 
for the equilibrium velocity v* in Eq. (11) needs to be 
modified. The parameters of order 0, 1, and 2 in the 
expansion are ho, hi, and ha. The definition of the ve- 
locities v*{£,t) and v(€,t) differ with the method used. 
We present four possible implementations which all as- 
sume /iQ = 0, since otherwise a source term in the mass 
balance would have to be taken into account. The sums 
in this paragraph run from i = 1 , . . . , A^ and the quan- 
tities hi, ha, fi, V, V*, and b arc functions of £ and t 
unless specified otherwise. 

The first method to implement a body force is referred 
to as METHOD A in the remainder of the paper. It uses 

^^=(^"^)^' h2=(^l-^)(v*b+bv*), (18) 

and a modified definition of v* and v which causes the 
infiuence of temporal and spatial derivatives of b on 
the density and momentum changes to vanish. For this 
method one obtains v* = v, with 



{Y.fiCi/j2fi)+Ath/2 



(19) 



instead of Eq. (5). A multiscalc expansion in time of 
the resulting discrete LB equation yields that the macro- 
scopic density p and velocity v recover the Navier-Stokes 
equations with an external body force term [35]. The 
forcing is applied in two steps during every time step 
At, one half within the collision step by the definition 
of V* and the second half within the streaming step by 
the term ipi. In the case of LB-MRT the part which is 
applied during the collision step is added to the modes 



mi{£,t) with i — 4,6,8, which represent the momentum 
flux. 

The second method (method B) is defined by setting 



hi = b, h2 = 0, 



(20) 



so that (pi{£-) does not depend on v*{£, t). The fufl accel- 
eration is apphed only within the streaming step through 
the term fi{t)- One sets 



= 5I/'c7E-/'*' 



(21) 



and the macroscopic velocity v defined as in Eq. (19). 
This simplification is useful because it reduces the com- 
putational effort, but it is restricted to stationary flows. 
In our simulations b(£) is time independent and we 
are mainly interested in the permeability and stationary 
flows so that we have adopted method B in our simula- 
tions below. In method B the macroscopic fields fulfill 
mass balance, but some additional unphysical terms ap- 
pear in the momentum balance [35]. Here we assume 
that all these additional terms are negligible or vanish 
for stationary flows, because we expect that all spatial 
gradients arc sufficiently small. 

METHOD C is intended for constant b and uses the 
same parameters hi and h2 as method B [36]. How- 
ever, the macroscopic velocity v = v* is calculated as 
in Eq. (21). This recovers momentum balance, because 
unphysical terms either vanish or are negligible, but it 
does not recover mass balance, which in this case reads 

|^+V.(pv) = -^V.(pb). (22) 

The reason is an inaccurate calculation of the macro- 
scopic velocity v(.£,i) [35]. The impact of this issue on 
the simulation results is shown in Sec. V. 

METHOD D suggests to incorporate the acceleration 
not by using the forcing term, but by adding the term 
Tb(£,t) to the equilibrium velocity v*(.£,t). The macro- 
scopic velocity v(.£,i) remains calculated by Eq. (5) [37]. 
This is equivalent to using the forcing term with 



hi 



rbb + bv* +v*b. 



(23) 



and V = V* given by Eq. (21). This implementation leads 
to the same drawback in the mass balance equation as in 

METHOD C. 

The most common boundary conditions (BC) used 
jointly within LB implementations are periodic (PBC) 
and no-slip BC. When using PBC, fluid that leaves the 
domain, i.e., the term l + /S.li in Eq. (8) exceeds the com- 
putational domain size, enters the domain from the other 
side. The no-slip BC, also called simple bounce-back rule 
(SBB), approximates vanishing velocities at solid sur- 
faces [13]. If the lattice point £-(-A£j in Eq. (8) represents 
a solid node, the discrete LB equation is rewritten as 



where the probability function /* is associated with c*, 
where c* = — c, is the probability function in opposite 
direction to /j. Midplane BC [37] improve the SBB 
eliminating the zig-zag profile when plotting the mass 
fiow g vs. ^3, but yield the same mass fiow Q, see 
Eqs. (27) and (28) for their definition, respectively. The 
SBB scheme depends on viscosity and relaxation time r. 
especially in under-relaxed simulations (large values of 
r) [36] . The numerically exact position of the fluid-solid 
interface changes slightly for different r which can pose a 
severe problem when simulating flow within porous me- 
dia, where some channels might only be a few lattice 
units wide. The permeability k, being a material con- 
stant of the porous medium alone, becomes dependent 
on the fluid viscosity. As demonstrated below within 
the LB-MRT method this k,-t correlation is significantly 
smaller than within LB-BGK [27, 38]. Recently, further 
improvements for no-slip BC have been discussed [33]. 
Most of these implementations use a spatial interpola- 
tion. For example, linearly and quadratic interpolated 
bounce-back [39, 40], or multirefiection [41]. To calcu- 
late boundary effects these methods use multiple nodes in 
the vicinity of the surface. For this reason these schemes 
are unsuitable in porous media where some pore throats 
might be represented by 2 or 3 nodes only. Consequently, 
we use midplane BC as well as PBC for our simulations. 
To drive the fiow on-site pressure or flux BC [42, 43] 
may be used. Using them it is possible to exactly set 
the ideal gas pressure (or density, see Eq. (6)) or flux on 
a speciflc node. Thus, creating a pressure gradient by 
flxing cither the pressure or the mass flux at the inlet 
and outlet nodes are feasible alternatives. 



III. SIMULATION SETUP 

The computational domain (see Fig. 1) £ is composed 
of three zones: the sample S describing the geometry and 
two chambers X (inlet) and O (outlet), before and after 
the sample, containing fluid reservoirs. The notation 



C{a) := {£ e £ : ^3 = a} 



(25) 



/*(£,t + At)-/i(^,t)=Atr!i(£,t), 



(24) 



denotes a cross-section, where C{Lx), C{Lx + l), C{L^ — 
Lq), and C{L'i—Lo + \) represent the cross-sections right 
before the sample (C(ii) S X), the first {C{Li+\) G S), 
and last (C(L3— Lq) € iS) cross-section within the sample, 
and the cross-section right after the sample (C(L3— Lq + 
1) € O), respectively. Every lattice point (node) in £ 
is either part of the matrix, denoted M, or part of the 
fluid, denoted V, so that MUV = C and MnV ^9. 
Results are presented in the dimensionless quantities 

x^x/Aa;, t = i/At, p = p/p\ p^p/{?,c,^p°), 

V = vAt/Ax, f = r/Af, fbuik = Tbuik/At, 

K = K/(Aa;)2, 6 = 6(Af)VAa;, q = qM/{p° {Axf), 



Computational domain C 



O 



^ 



H 



K 




Lt 



FIG. 1: The computational domain £. The (porous) sample 
is 5, and the fluid is accelerated in the acceleration zone T. 
Two fluid chambers T and O are used to avoid artifacts. 



where the discretization parameters Ax and Ai are cho- 
sen according to the analysis presented in Sec. II. Unless 
otherwise noted, the relaxation time is f = 0.857 and 
for LB-MRT simulations fbuik = 1-0 is used. Generally, 
results from LB simulations are labeled with the super- 



script 



e.g., the density p""^. If the results refer to a 



specific implementation (BGK or MRT) they are labeled 
accordingly, e.g., p^^^ or p^^"^ . 

The fluid is driven using MODEL B. The acceleration 
b = 6 es is not applied throughout the whole domain but 
only within the acceleration zone J- <Zl. An acceleration 
oih = 10~^ is used for all simulations. 

The average for a physical quantity w is 



mE-W' 



IV 



(26) 



t.ev 



with the domain V £ {C,S,V,I,0,F,C{a)} and |V| 
the number of nodes in that domain. The mass flow 
q through a cross-section C{a) is given by 



q{a)= J2 pW^3W(Ax)2 



(27) 



eec{a)nv 



with p[€)Vy,{(.) being the momentum component in di- 
rection of the flow. The mass flow through the whole 
domain is 



1 ^' 



IV. CALIBRATION 



(28) 



To calibrate the simulation we simulate Poiseuille flow 
in pipes with quadratic cross-section. The simulation 



parameters are defined by 



2 
L2 + I 



Axei 



Aa;e2+ [i^--] Axes, 



S = {£e/::2<£i <Li-l, 

2 < ^2 < ^2 - 1, 4 < ^3 < L3 - Lo}, 



O = {£eC 
T = {£eC 



h < Li}, 

L3-Lo<e3<L3}, 
^3<2} 



(29) 



where Lj = 3, Lg = L3 — 6 and Lo = 3. The system 
dimensions are Li ~ L2 = B+2, L3 = AB, with B/Ax = 
B e {4, 8, 16, 32, 64} the channel width. 

According to Rcf. [44] the analytical solution for the 
velocity component in flow direction in a pipe with 
quadratic cross-section is 



w™(xi,.T2)= lim v{xi,X2,M) 



(30) 



,/(2n+l)7r \ /(2n + l)^ 

cosh I :; Xi COS X2 



Cn = (-1)"- 



B 



B 



{2n + 1)3 cosh 



(2ri + l)7r 



where xi £ [-B/2,B/2],X2 G [-B/2,B/2]. The Carte- 
sian coordinates xi and X2 have their origin in the cen- 
ter of the pipe, {^p)^ is the pressure gradient in flow 
direction and 77 the dynamic viscosity. The expression 
v{xi,X2,M) is asymmetric in xi and X2- Contrary to 
the no-slip condition the velocities v{B/2, X2,M) are not 
zero for finite M. To estimate the truncation error we 
define 



v{x,M) 



277 



(yp),B^ 



v{B/2, X, M), X e [-B/2, B/2] 



and 



(31) 



K'wall(M)||. 




?/2 



\vix,M)f dx, (32) 



-B/2 



with v(x, M) being the normalized velocities on the wall 
calculated from Eq. (31). ||'ywau(-^)|l2 quantifies the 
truncation error at finite M . Requiring that the trun- 
cation error ||wwaii(-^^)|l2 ^^ ^^ least three to four decades 
smaller than the velocities in the corners, for example 
w(g(l, l,i3/2)i,g(l, 1,^3/2)2, Af) or any other such cor- 
ner velocity, yields M « 50. For all further comparisons 
with LB simulations we use M = 200. If M is chosen too 
small, a meaningful comparison of the simulation results 
with the analytical solution is not possible because of the 
inaccuracies in the numerical evaluation of the analytical 
solution itself. 



Eq. (30) is the stationary solution for the velocity com- 
ponent in flow direction on a quadratic cross-section in 
an infinitely long pipe and for a constant pressure gra- 
dient. Therefore, the simulated w^<^^(£,i) and p=^°'^(^,t) 
are inspected for convergence at the end of the simula- 
tion t = fend and the assumption of a constant pressure 
gradient is checked. We define 



Sw(t, dt) = max 
t.£(snv) 



w{£,t)-w{jE,t~dt) 
w{i,t) 



(33) 



^'="'(£,f) w 10"'' so the absolute 
^^, using the relative changes Sv 

^BGK/ 



as the maximum relative change of a quantity w dur- 
ing the time dt and within the computational domain 
iS n "P, where w{i,t) is either the velocity v^'^'^{£,t) or 
the density p^°^{(.,t). Because the pressure is propor- 
tional to the density, Eq. (6), the pressure is converged, 
if the density is sufficiently converged. The results from 
Eq. (33) are shown in Tab. I. In the simulations the ve- 
locities are of order v 
changes are of order 10 

from Tab. I. The fiuid density is p^'''^{e.,t) w 1.0 giving 
absolute changes of order 10~^ . The variation of the pres- 
sure gradient can be approximated by 2Sp/{L3) < 10~^. 
When calculating errors by comparing them with ana- 
lytical solutions the number of significant digits is de- 
termined by the convergence of the simulation. We use 
the notation u^<='«(£), p^'='^(£) and p^<='«(^) for the ve- 
locity, density and pressure at the end of the simulation 

f — fend- 

Due to the way we drive the fiow, the pressure in- 
creases in the acceleration zone J-' and then decreases 
along the fiow direction. See Fig. 8, where the average 
density (/5)g(^, ^pp — 1, from a LB-BGK simulation for a 

pipe of width _B == 7 is shown. To verify that the pres- 
sure gradient (Vp)^ can be assumed to be constant as 
required by Eq. (30), we linearly fit {p^'^'^)c{e3)nv inside 
the sample. In the LB simulations, for pipes of width 
B = 4, 8, 16, 32, and 64, all residues of the linear fit are 
of the order 1 x 10~^, so that the pressure gradient can 
be assumed to be constant. 

Next, the velocity component in fiow direction v^^{£) 
with £ £ C{L^/2) is compared to the analytical 
solution Eq. (30), evaluated at the node positions 
w'^"(g(£)i,g(£)2). The cross-section C{L'i/2) is chosen 



B 


dt/At 
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2.34 


0.174 
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1 


20 


2.80 


0.174 


16 


5 


30 


4.27 


0.869 


32 


5 


50 


7.45 


0.869 


64 


10 


120 


1.33 


1.74 



TABLE I: Maximum relative change of the velocity 
Sv{tend, dt) and density 5p(iend, df), Eq. (33), during the time 
dt when the simulation ended at tend. B is the dimensionless 
channel width. 



to minimize finite size effects and artifacts from the 
in/outlet chamber. We define absolute and relative errors 
of the velocities as 
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Fig. 2 provides an overview on the structure of |e^'^'^(£)| 
and Fig. 3 shows |e^'^'^(£)| with £2 = L2/2 and £3 = L3/2 
as a log-linear plot for different pipes of width B and 
LB-BGK. The largest relative errors are located in the 
corners and close to the wall. As the resolution increases 
the relative error declines rapidly. In the central region 
it is much smaller than 1%. To gain insight into how 
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FIG. 2: Overview of the relative error \e^°'^{£)\ with £3 = 
L3/2, for pipes of widths B G {4, 8, 16, 32, 64}. Nodes at the 
corners cause the largest error followed by those close to the 
solid walls. For larger pipes the error decreases substantially. 



strong the value of the relaxation-time r influences the 
accuracy of the velocity field, simulations with LB-BGK 
and LB-MRT and different relaxation-times f = 0.7, 
1.0, 2.0, 2.5, and 3.0 are investigated. Fig. 4 displays 
the relative error of the velocity e]f{£), Eq. (35), with 
I2 = B/2, £3 = L3/2, and B = 20 for both implemen- 
tations LB-BGK and LB-MRT and different relaxation 
times. The calculated velocity in the center of the pipe 
is in good agreement with the theoretical solution, hav- 
ing a relative error smaller than 1%. It is interesting to 
note that when using the LB-BGK method the largest 
error occurs for a large relaxation time f = 3.0 (over 
relaxation), whereas the largest error for the LB-MRT 
result occurs at a small relaxation time f = 0.7 (under 
relaxation). The calculated velocities tend to be overes- 
timated ey'^^{£) > for LB-BGK simulations and under- 
estimated ej;™'^(.£) < for LB-MRT simulations. When 
using LB-MRT the relative error is smaller by roughly 
a factor 10~^ when compared to results of the LB-BGK 
method. 
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FIG. 3: Log-linear plot of the relative error along the central 
line on a cross-section, i.e |e^'^'^(£)| with £2 = L2/2 and I3 — 
I/3/2. The error is largest at the walls and declines towards 
the center of the pipe. Different line styles indicate different 
pipe widths B = 4, 8, 16, 32, 64 as shown in the legend. 



For permeability calculations from Darcy's law, see 
Eq. (36), the mean velocity {v^^)s is used. Therefore, 
the mean relative error (jeJ^'^'^Ds and the mean absolute 
error (|e^°'^|)5 are of interest. Both decrease when B 
increases, as shown in Tab. IL The mean relative er- 
ror shows a power law behavior (|e^'^^|)5 = oiB"^, with 
parameters ai ~ 0.6 and 02 ~ —1.6. This relation can 
be used to calculate the relative error for arbitrary pipe 
widths B. Overall, the LB-BGK implementation is able 
to reproduce the velocity field for quadratic pipes very 
accurately. The mean relative error (jeJ^'^'^Ds is below 
1% if the pipes are resolved better than B > 14. 

Following the evaluation of the calculated velocity 
field, permeabilities are calculated using both implemen- 
tations LB-BGK and LB-MRT. The permeability k^'^ is 
calculated using Darcy's law: 



MRT 



-v 



((Vp) 



LB\ ■ 

3 /S 



(36) 



where 77 is the dynamic viscosity, {v^^)g is the average 
velocity in the sample and {(Vp)^ ) c is the average pres- 
sure gradient component in direction of the flow. Details 
how {{Wp)^^)g can be determined from p^^{£) will be 
discussed later in this article. The dynamic viscosity is 
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(|er"l).,[xl0-«] 
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0.064099 


0.300 
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0.021840 


0.114 


16 


0.007241 


0.068 


32 


0.002331 


0.034 


64 


0.000784 


0.023 



TABLE IL Mean relative error (|e?'^'^|)5 and absolute error 
(Ifiu'^'^Ds for different pipes with width B. The error declines 
rapidly as the pipe width increases. 
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FIG. 4: Relative error e^^(.«) with ^2 = B/2 and £3 = L3/2 
for flow in a quadratic pipe of width B = 20 and different 
values of f . The upper figure shows the results for the BGK 
implementation and in the lower figure the MRT results are 
shown. The relative error of the MRT simulations is smaller 
by 10"^ than the BGK results. 



calculated as 7] — vp^^ with p^ 
mated by 



'\P^'')vns approxi- 
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The analytically obtained permeability is [44] 

B^ (1 64 ^tanl^ ((2^^+1)1)' 
K™(m= lim — --^V- ^ 2/ 



(2n-fl)5 



(37) 



(38) 



where we use Af = 200 for numerical evaluation. To 
evaluate the error we define 






K™{B) 



(39) 



The relative errors e^'^^{B) and e^^'^{B) are shown in 
Fig. 5 and it can be observed that they fall below 1% 
for all pipes wider than B = 16Aa;. It seems that the 
LB-BGK method is slightly more accurate, but the re- 
laxation time f = 0.857 was fine tuned to reproduce the 



exact result with LB-BGK. Fig. 5 shows that an adjusted 
relaxation parameter f can make up for the methodically 
inferior LB-BGK implementation. In realistic porous me- 
dia, however it is not possible to determine an optimal 
relaxation time r, because the pore diameters and pore 
throats vary, although a useful range of r can be deter- 
mined, see Sec. V. Therefore the LB-MRT method is 
more reliable as its results are less dependent on r. 

The results for the velocity field and permeability show 
that even for a simple quadratic channel a resolution of at 
least 20 lattice nodes is required to achieve an accuracy of 
the permeability of order 1%. At present, discretization 
at this resolution is neither experimentally available nor 
computationally manageable. 




FIG. 5: Relative error €^,{3), of the permeability k^^{B) ver- 
sus channel width B at fixed resolution Ax — 10~^m as calcu- 
lated using LB-BGK (solid line) and LB-MRT (dashed line) 
simulations. 
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FIG. 6: Relative error ej:;^ vs. the value of f for a Poiseuille 
flow in a quadratic pipe with different pipe width B = 
1, 2, 3, 4, 5, 10 as indicated by different line styles. 
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V. POTENTIAL DIFFICULTIES LEADING TO 
INACCURACIES 

In this section we discuss typical difficulties arising 
when calculating permeabilities for complex geometries. 
This includes the influence of the relaxation time r on 
the permeability, the accurate approximation of the aver- 
age pressure gradient, the implementation of the external 
force and the discretization error. 

When using SBB, the relaxation time r slightly 
changes the position of the boundary between adjoined 
fluid-solid nodes. Due to this effect the relaxation time 
has a substantial influence on the permeability calcula- 
tion [33]. One also has to be aware that this effect is 
always correlated with the discretization error and can- 
not be corrected analytically when investigating stochas- 
tic porous media. To analyze the influence of r on 
the permeability we again investigate Poiseuille flow in 
quadratic pipes using a computational grid aligned with 
the pipe geometry to minimize the discretization error. 
The relative error e)^^(T) is shown in Figs. 6 and 7. For 
over-relaxed systems (f > 1), the LB-BGK method yields 



FIG. 7: Relative error ^ for Poiseuille flow in a quadratic 
pipe with width B = 10, calculated using LB-BGK, dashed 
line, and LB-MRT, solid line. On the left, the interval f £ 
[0.5,4.0] is shown and on the right, a zoom in the interval 
f £ [0.5, 1.0] can be seen. 



incorrect results, increasing dramatically the dependence 
of permeability on r when the geometry is poorly dis- 
cretized {B < 5). If f € [0.5, 1.0], the absolute error of 
permeability estimation is less than 3% for all B > 2. 
The LB-MRT method has to be considered more reli- 
able in general because the influence of r is much smaller 
and the influence of Tbuik is practically insigniflcant. For 
example, in Fig 7, using a value f = 3.5, the error of 
the LB-BGK method is 53.70%, while the error of the 
LB-MRT method is 4.213%^ However, if only a smaU 
number of nodes is used {B > 5) even the LB-MRT 
method produces a substantial error. Fig. 7 compares 
the results between the LB-BGK and LB-MRT for a 
pipe with B = 10. Here, both absolute errors in the 
interval f e [0.5, 1.0] are smaller than 1.5%. It is impor- 
tant to stress here that outside the interval f S [0.5, 1.0] 
the LB-MRT results remain accurate when B decreases. 



which is not the case for LB-BGK. 

To compute the permeabihty using Darcy's law as 
given in Eq. (36), the average pressure gradient in di- 
rection of the flow ((Vp)3)g is required. Because the 
permeability can strongly depend on the way the pres- 
sure gradient is obtained, alternative methods for its de- 
termination are discussed: 

a) Calculating the slope of a linear fit through the full 
data set {p)^/^ j^p obtained using all cross-sections 

C(£3), iI + l<^3<i3-io• 
b) As a), but using only the cross-sections C{i,z), Lx + 
1 + W <i3<L3-Lo-W,W elN,see Fig. 8. The 
cross-sections closer than W to the inlet and outlet of 
the sample are not taken into account. The idea is to 
minimize boundary effects. 

c) Approximation of {{^p)^)^ by the arithmetic mean of 
the pressure at C{Lx + l) and C{L3 — Lo), i.e., 
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3/S 



{Ls - l)Ax 
where Lg is the sample length. 



(40) 



d) Approximation of {i^p)^)^ by the arithmetic mean of 
the pressure at C{Lx) and C[L'i, — Lq -1-1), i.e.. 
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C(L3-Lo + l)r\V 
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c(Lx)nv 



3'S 



{Ls + l)Ax 



(41) 



For a quantitative comparison of the different methods 
simulations with the following parameters are performed: 

ii = 40, L2 = 40, L3 = 80, 

S = {£ e C : 17 < ei < 23,17 < £1 < 23, 21 < £3 < 60}, 
Z = {£ e £ : 1 < £3 < 20}, O = {£ e £ : 61 < £3 < 80}, 
J" = {£ e C : 6 < £3 < 15}, Lx = 20, L5 = 40, Lo = 20. 

The average density (p)c(£ -jp-p, is shown in Fig. 8. Al- 
though the density field p{i) is continuous, the average 
ip)c(t )nv shows two discontinuities, one at the beginning 
of the sample {£3 = 20) and one at the end of the sample 
{£3 = 60). These can be explained by the small compress- 
ibility of the fluid. The majority of the fluid in chamber 
I flows towards the surface of the sample causing an in- 
creased local density. The same effect can be observed 
right behind the sample where one finds a low density 
due to the fluid compressibility. Because we use periodic 
boundary conditions, the pressure is almost constant in 
both chambers X and O and only increases in the accel- 
eration zone J-'. The small increment right before the 
sample and the small decrement right after the sample 
are both imperceptible in Fig. 8. The results presented in 
Tab. Ill show that by using alternatives a), b) and c) an 
error smaller than 1% can be obtained. The remaining 
method d), however, shows a substantially larger error 
and is therefore not suitable for measuring the pressure 



gradient. Alternative b) is not taking into account the 
cross-sections closer than 1^ = 10 to the inlet and outlet 
of the sample. Changing W does not influence the ac- 
curacy much. For stochastic porous media we suggest to 
use alternative c) because it is very easy to implement 
and no fit is necessary. The last row of Tab. Ill shows 
the values obtained without an injection chamber and 
with a force acting throughout the whole domain. Even 
though the result is accurate, this method has a major 
disadvantage, because it can only be applied to periodic 
samples and not to stochastic porous media. In realistic 
porous media, chambers before and after the sample are 
necessary to provide a fluid reservoir but they might de- 
crease the accuracy of the method due to disturbances of 
the velocity fleld at the in- and outlet. 



B = 7 


a) 


b) 


c) 


d) 


k'-'' 


1.73244 


1.73689 


1.71346 


1.64270 


e. [%] 


0.61423 


0.87278 


0.48787 


4.59761 


Without injection channel 


ft^^ 


1.73689 


e. [%] 


0.87277 



TABLE III: Results of the three alternatives to measure 
((Vp)3)„ and the case without using an injection channel. 
Shown are calculated permeabilities and their relative error 
for B = 7. For alternative b) H^ = 10 is used. 




as obtained 



8: Average density {(p)c(e^)nv ~ 1) v' 

The sample S is placed in 



FIG. 

from a LB-BGK simulation. 
£3 £ [21,60], leaving 20 nodes before and after the sample 
as injection chamber. Within the acceleration zone J-" at 
£3 £ [6, 15] the density and pressure increase. Alternatives 
to measure the density or pressure gradient: a) Use all nodes 
inside the sample for a linear fit (0). b) Linear fit not using 
the nodes closer than ly = 10 to the point where the fluid 
enters or leaves the sample (x). c) Use the first and last 
cross-section C{Lx+l) and C^L^—Lo) inside the sample (+). 
d) Use the cross-section C{Lx) and C(I/3 — Lo -I- 1) (■). The 
nodes represented by D define the injection channel (X and 
O). 

An important point for performing high precision per- 
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mcability measurements is the way the pressure gradi- 
ent is generated. While pressure boundary conditions 
provide a well defined way of fixing the pressure at the 
in- and outlet, they assume an ideal gas and are slightly 
harder to implement than a simple body force driving the 
flow. In addition, even though the pressure is fixed be- 
fore and after the sample, an injection chamber is still re- 
quired and for high precision permeability measurements 
one has to measure the pressure gradient as discussed 
above. Therefore, most LB implementations found in 
the literature use body forces. In fact, all papers we are 
aware of, that have been published before 2002, and a 
large fraction of more recent publications use an incor- 
rect force implementation which can lead to severely er- 
roneous permeabilities. Popular examples for such imple- 
mentations are method C and method D. They lead 
to an underestimation of the velocity v(£) in the direction 
of the flow on the lattice nodes where the acceleration is 
acting. Many publications apply the force throughout 
the whole simulation domain. The results obtained from 
such implementations cannot be trusted for two reasons: 
Firstly, in METHOD C and METHOD D the macroscopic 
velocity in the acceleration zone is smaller than the cor- 
rect value. In some cases it can even be negative. Sec- 
ondly, the pore structure plays an important role. The 
number of nodes at any cross-section €{£3) determines 
the number of times the additional acceleration term has 
to be added to the mass flow in order to assure a constant 
flux. 

In Figure 9 we compare method B and method C. 
All simulation parameters except the pipe width (B ~ 5) 
are kept as before so that 

S = {£eC:18< ii< 22,18<^2< 22,21< ^3< 60}. (42) 

The line representing the application of method C and 
the external acceleration applied throughout the whole 
domain has discontinuities exactly at the position where 
the width of the channel changes abruptly, i.e. at £3 = 20 
and ^3 = 60 (local porosity dependency). The line rep- 
resenting the application of method C throughout an 
acceleration zone !F shows discontinuities within the ac- 
celeration zone, i.e. in the interval €3 € [6,15]. These 
discontinuities are not present when using method B. 
The term rj{v^)s in Darcy's law (see Eq. (36)) is usually 
approximated by Qv / A, where the total mass flux Q is 
calculated by averaging 9(^3) in the whole sample. A 
represents the sample cross-sectional area. If an external 
acceleration is not implemented correctly, the calculated 
Q is always incorrect leading to a wrong estimate of k. 
For the example in Fig. 9 which uses method C and a 
force throughout the whole domain, Q is underestimated 
but remains positive. This is not the case if an accelera- 
tion zone F is used. Here, the calculation of Q leads to 
an unphysical negative result, so that the permeability is 
always underestimated and in some cases negative. Such 
cases can also be observed for inhomogeneous stochastic 
porous media, where the variation of pore sizes is very 
large [45]. 
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FIG. 9: Mass flux qil-j) in a quadratic pipe of width B = 5 
as obtained from a LB-BGK simulation. If the acceleration is 
implemented as deflned by method C, 9(^3) is not constant in 
the regions where the acceleration is applied, method B en- 
sures a correct constant mass flux 5(^3) throughout the whole 
domain. 



Another important issue is the effect of discretization. 
When investigating square pipes the lattice is aligned 
with the solid-fluid interface. This is not the case for 
the simulation of flow in realistic stochastic porous me- 
dia. Thus, the influence of discretization effects is sub- 
stantially larger than in the ideal cases presented before. 
We investigate the order of the resulting error by calcu- 
lating the permeabilities in pipes with a circular and an 
equilateral triangular cross-section. The samples are of 
size Li ~ J, L2 = J and. L3 = 4J with J € N and 
the cross-sections are defined by their diameter Bq (cir- 
cular) or their side-length B/x (equilateral triangle) with 
Bq — _Ba = (>/ ~ 2)A.T. The analytical solutions for the 
permeabilities of circular and equilateral triangular pipes 
are 
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(43) 

(44) 

Dis- 
cretizing these areas on a cubic lattice results in approx- 
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where A™ and A™ are the cross-sectional areas. 



imate cross-sectional areas A'^ and A^. 



Let 



and 



^o 



be the relative discretization errors of those areas. 
The permeabilities as calculated from the simulation re- 
sults are k^^ and n^^, with their relative errors be- 



ing 



,BGK 



and 



e.. 



Fig. 10 depicts that for both ge- 
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omctries the relative error of the permeabilities is much 
larger than for pipes with quadratic cross-section, see 
Fig. 5 for comparison. Furthermore, it can be seen that 
the error in permeability correlates with the error of the 
discretized area. This discretization error is not present 
when investigating square pipes that are aligned with 
the grid. In stochastic porous media this discretization 
error is inevitable. Therefore, arbitrarily structured pore 
throats have to be resolved at a much higher resolution 
for high precision permeability calculations. This is a se- 
rious limitation when calculating permeabilities for labo- 
ratory sized porous media using the techniques discussed 
in this article. 




FIG. 10; Relative errors in the area discretization t^^ and 
permeability estimation ej?'^*' for a circular and triangular 
cross-section pipe with different system sizes J. The per- 
meability error correlates with the discretization error. Com- 
pared to the results for quadratic pipes, the permeability error 
is at the same resolution approximately twice as large. 



VI. APPLICATION TO FONTAINEBLEAU 

SANDSTONES 



After validating the simulation results, determining er- 
rors for pipe flow, and pointing out problems when calcu- 
lating permeabilities with LB implementations, we now 
apply our findings to investigate a porous sample. We 
calculate the permeability of a sample of Fontaincblcau 
sandstone, gained by thresholding a discretized /x-CT 
data set. This particular data is chosen, because it 
has been investigated previously using a finite difference 
method and another LB-implementation [12]. The re- 
sults for the calculated permeabilities in [12] show ex- 
cellent agreement with the experimental results in [46], 
therefore this sample serves as a benchmark for the per- 
meability calculations presented here. Calculations are 
carried out at a low (l.r.) and high resolution (h.r.). The 



l.r. computational domain is 

Li = 305, L2 = 305, ig = 320, 

5 = {£e£:3<£i< 303, i < ^2 < 303, 11 < £3 < 310}, 
I = {£ e £ : 1 < £3 < 10}, O = {£ e £ : 311 < 4 < 320}, 
J" = {£ e £ : 1 < 4 < 5}, ii = 10, L5 == 300, Lo = 10, 

Aa; = 7.5x10"%!. 

The high resolution sample is created from the low resolu- 
tion sample by substituting every voxel with eight voxels 
on a cubic sublattice. The h.r. computational domain is 

Li = 605, L2 = 605, L3 = 620, 

5 = {£e/::3<^i< 602, 3 < ^2 < 602, 11 < £3 < 610}, 
X = {£ e £ : 1 < ^3 < 10}, O = {£ e £ : 611 < ^3 < 620}, 
J" = {£ G £ : 1 < ^3 < 5}, ii ^ 10, L5 = 600, Lq = 10, 

Ax = 3.75xl0~*^m. 

For the permeability calculation an approximation of 
Darcy's law is used, see Eqs. (36), (40) and (37), 



i/(L5-l)Aa;(w3)5((p)c(L^+i)nP + (/5>c(L3 
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(45) 

The kinematic viscosity v is calculated using Eq. (15) 
with f, as in Tab. IV, fbuik = 1-0 and At from Eq. (7) 
with Cg ~ lms~^. Simulations arc performed with the 
LB-BGK and LB-MRT method for 100 000 time steps. 
The relaxation times used in the simulation and the 
calculated permeabilities for the l.r., h.r. sample and 
from [12] are given in Tab. IV. The calculated permeabil- 
ities Ki.r. and Kh.r. were linearly extrapolated for infinite 
resolution at 1/ Lg ~ yielding Kextrap- 

Our resuh k^.,. = 608 [mD] for f = 0.688, see Tab. IV, 
is in good agreement with the result in [12] being k ~ 
621 [mD]. Although the simulation setup is different, a 
relative difference of only 2% is obtained. Fig. 11 (top) 
confirms that permeability results gained from LB-BGK 
simulations are particularly dependent on the relaxation 
time f that is used. Therefore, when investigating com- 
plex geometries, where f cannot be optimized for a spe- 
cific geometrical shape, a LB-MRT should be used. As 
expected, when calculating permeabilities for complex 
geometries, the influence of f is much stronger than 
within simple geometries, i.e., square pipes, see Sees. IV 
and V together with Fig. 7. The extrapolated permeabil- 
ities Koxtrap are an estimate for the true permeability of 
the discretized /i-CT sample at resolution 3.75 x 10~^m 
and not the true permeability of the sandstone. To esti- 
mate the true permeability of the sandstone by extrap- 
olation, new /i-CT data with higher resolutions would 
be required. However, from Fig. 11 (bottom) it can be 
seen that the extrapolated permeability values have a 
small spread, in a range from 473-510 [mD] regardless of 
the simulation method and relaxation time used. This 
indicates that, if sufficiently high resolved sample data 



12 



and computer performance arc available, an extrapola- 
tion analysis, even using LB-BGK f = 1 results, might 
give a good approximation of the true permeability. One 
could expect that errors due to random geometries would 
experience random cancellations. However, the results 
presented in Fig. 11 clearly show that this is not true for 
a realistic porous medium. 



if applied to the Fontainebleau sandstone (see Fig. 11). 



LB Method 


f 


Kl.r. [mD] 


Kh.r. [mD] 


«:extrap [luD] 


BGK [12] 


0.688 


621 


— 





BGK 


0.688 


608 


559 


510 


BGK 


0.857 


773 


634 


495 


MRT 


0.688 


505 


489 


473 


MRT 


0.857 


558 


518 


478 


MRT 


1.000 


601 


541 


481 



TABLE IV: Simulation results of the permeability calcula- 
tions for the Fontainebleau sandstone. Results for the low 
resolution sample are labeled Ki.r. , the high resolution sample 
results Kh.r. and the extrapolation results Kextiap. Different 
relaxation-times f and LB-BGK and LB-MRT implementa- 
tions are compared. 



VII. CONCLUSION 

Our simulation setup of an acceleration zone and 
in/outlet chambers, together with our approximations of 
Darcy's law, provides a method for permeability calcula- 
tions. Several problems in the numerical implementation 
and data evaluation were addressed, such as a correct 
acceleration implementation and an adequate approxi- 
mation for calculating the pressure gradient. Caveats 
when using LB simulations to calculate permeabilities 
have been exposed. Wc performed detailed studies with 
different LB-implementations, i.e., BGK and MRT, and 
for various systems to quantitatively determine the accu- 
racy of the calculated velocity field and calculated per- 
meability. We find that for reasonably resolved quadratic 
pipes, the error of the calculated permeability is be- 
low 1%. Investigating non-aligned geometries, circular 
and triangular pipes, the discretization and permeabil- 
ity error is roughly 4% at comparable resolutions. From 
this we infer that permeability calculations in stochastic 
porous media will have a significantly larger error, be- 
cause the resolution of pores and pore walls is usually well 
below the resolution used for our pipe calculations above. 
Comparing the two LB-implementations, LB-BGK and 
LB-MRT, we find that LB-MRT reduces the dependence 
of the permeability on the value of t substantially. Using 
LB-BGK and a relaxation time t tailored to give good 
results for a specific geometry docs not assure reliable 
results in a stochastic porous medium. For example, we 
found that LB-BGK and f = 0.857 yields the best result 
for 3D Poiseuille flow in a quadratic pipe. However, for 
this value LB-BGK and LB-MRT resuhs differ by 20% 
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FIG. 11: The top plot shows the permeability for different f, 
LB schemes, the l.r. and the h.r. sample. The influence of f on 
the permeability is stronger using the BGK implementation. 
The bottom plot shows the permeability results, for the l.r. 
sample with Ls = 300 and the h.r. sample with Ls = 600. At 
IjLs = the extrapolated permeabilities ftextrap are shown. 



Therefore, LB-MRT is suggested to be used for perme- 
ability estimates based on LB simulations. Further inves- 
tigations using the LB method for flow through stochas- 
tic porous media should include a resolution and relax- 
ation time dependent analysis together with an appropri- 
ate extrapolation scheme for more reliable permeability 
estimates. 
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